library(tidyverse)
library(DeclareDesign)
library(rdss)

sims <- 2000
bootstrap_sims <- 2000

N <- 1000
ate <- 0.20

design_1 <-
  declare_model(
    N = N,
    U = rnorm(N),
    potential_outcomes(Y ~ ate*Z + U)
  ) +
  declare_assignment(Z = complete_ra(N)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z)

studies_1 <- simulate_design(design_1, sims = sims)

studies_1_for_merge <- 
  studies_1 |>
  transmute(sim_ID, 
            original_estimate = estimate,
            original_std.error = std.error,
            original_conf.low = conf.low,
            original_conf.high = conf.high)

designs_2 <- redesign(design_1, N = seq(2000, 10000, 1000), ate = 0.15)

studies_2 <- simulate_designs(designs_2, sims = sims)

is_equivalent <-
  function(estimate, std.error, tolerance) {
    p_upper_eq <-
      pnorm(q = (estimate + tolerance) / std.error,
            lower.tail = FALSE)
    p_lower_eq <-
      pnorm(q = (estimate - tolerance) / std.error,
            lower.tail = TRUE)
    (p_upper_eq <= 0.05) & (p_lower_eq <= 0.05)
  }

simulations <- 
  studies_2 |>
  left_join(studies_1_for_merge, by = "sim_ID") |> 
  mutate(
    replication_outside = !(original_conf.low < estimate & estimate < original_conf.high),
    original_outide = !(conf.low < original_estimate & original_estimate < conf.high),
    abs_diff = abs(original_estimate - estimate),
    se_diff =  sqrt(original_std.error^2 + std.error^2),
    t_stat = abs_diff / se_diff,
    dis_sig = pnorm(t_stat, lower.tail = FALSE) * 2 < 0.05,
    dis_nonequiv = !is_equivalent(estimate = abs_diff, std.error = se_diff, tolerance = 0.2)
  )

write_rds(simulations, file = "diagnosis_objects/diagnosis_23a.rds")
